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We present measurements of the specific heat coefficient 7(= C/T) in the low temperature limit as 
a function of an applied magnetic field for the Fe-based superconductor BaFe2(Aso.7Po.3)2. We find 
both a linear regime at higher fields and a limiting square root H behavior at very low fields. The 
crossover from a Volovik-like i/H to a linear field dependence can be understood from a multiband 
calculation in the quasiclassical approximation assuming gaps with different momentum dependence 
on the hole- and electron-like Fermi surface sheets. 



INTRODUCTION 



The symmetry and detailed structure of the gap func- 
tion in the recently discovered iron pnictidei and chalco- 
gcnido 2 - high temperature superconductors is still under 
discussion. Across an increasingly numerous set of ma- 
terials families, as well as within each family where su- 
perconductivity can be tuned by doping or pressure, ex- 
perimental indications are that there is no universal gap 
structure!^ Instead, the superconducting gap appears 
to be remarkably sensitive to details of the normal state 
properties. This "intrinsic sensitivity"— may be due to 
the unusual Fermi surface topology, consisting of small 
hole and electron pockets, and to the probable A% g sym- 
metry of the superconducting gap which allows a contin- 
uous deformation of the order parameter structure from a 
fully gapped system to one with nodes (for a review see, 
e.g. Ref. 0). It is important to keep in mind, though, 
that another possibility to account for the observed vari- 
ability is that different experiments on the same material 
may probe selectively different Fermi surface regions and 
hence different gaps within the system. 

The Ba-122 family of materials has been intensively 
studied because large high quality single crystals are rel- 
atively easy to produce^i Within this family, the isova- 
lently substituted system BaFe2(Asi_ K Pa;)2 with a max- 
imum T c of 31 K is particularly intriguing because it 
exhibits a phase diagram and transport properties re- 
markably similar to the heterovalently doped system 
Ba(Fci_ x Co x )2As2 and displays many signatures of ap- 
parent quantum critical behavior at optimal dopingj^r— 
In the superconducting state, penetration depth^ NMR 
spin-lattice relaxation^ thermal conductivity temper- 
ature dependence^ and thermal conductivity angular 
field variation^ show clear indications of nodal behavior. 
Surprisingly, a linear field dependence of the specific heat 
Sommcrfcld coefficient 7 was measured^ on optimally 
doped samples from the same batch. Such a behavior is 



expected for a fully gapped single band superconductor 
since the fcrmionic excitations from the normal cores of 
vortices provide the only contribution to 7 at low T, and 
the number of these vortices scales linearly with the field 
H . It was argued in Ref. [TH that the specific heat mea- 
surement might be consistent with the other experiments 
suggesting nodes if the heavy hole sheets in the material 
were fully gapped, while the gaps on the lighter electron 
sheets were nodal. In such a case the 7 ~ \f~tl behavior 
would be difficult to observe in experiment. 

In this paper, we report new experimental data on the 
magnetic field dependence of the specific heat of opti- 
mally doped BaFe2(Asi_ 2; P a ;)2 samples. We have ex- 
tended our previous measurements to 15 T to higher 
fields up to 35 T (« |if C 2(0)), where we find a contin- 
uation of the linear behavior reported earlier. However, 
more precise measurements at low fields have revealed 
the presence of a Volovik-like s/H term which persists 
roughly over a range of 4 T, crossing over to a linear be- 
havior above this scaled The observation of this term, 
consistent with nodes in the superconducting gap, there- 
fore supports claims made in earlier workji2r— without 
the need to assume an extremely large mass on the hole 
pockets. 

Theoretical estimates using the Doppler shift method 
for isotropic gaps given in Ref. [TH were oversimplified, 
but did show the need for a more thorough analysis of 
anisotropic multiband systems, and stimulated further 
experimental work, both of which we report here. The 
theoretical difficulties can be seen easily by considering 
a simple two-band model with two distinct gaps Aj and 
A2; where we assume for the moment that A2 > Ai. If 
the two bands are uncoupled, the two gaps correspond 
to two independent coherence lengths & ~ v F,i/ '("Aj), 
where i — 1,2, and two independent "upper critical 
fields" -ff C 2,i- Vortex core states of the large gap A2 are 
confined to cores of radius ~ £2 ■ For fields in the range 
H C 2,i ^ H < if C 2,2, the vortex cores of the small gap 
will overlap, while the large gap cores will still be well 
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separated. Note that if Ai is very small (these consider- 
ations also describe crudely nodal gaps), this field range 
can be wide and extend to quite low fields. On the other 
hand, methods of studying quasiparticle properties in su- 
perconductors are typically adapted to calculating near 
H c i or H C 2, i.e. in the limit of widely separated or nearly 
overlapping vortices. The current problem apparently 
contains elements of both situations. In the absence of 
interband coupling, of course, one can use different meth- 
ods, corresponding to the appropriate field regimes, for 
the distinct bands. For coupled Fermi surfaces, however, 
such an approach is not viable. In the immediate vicinity 
of the transition, where the Ginzburg-Landau expansion 
is valid, there is a single length scale controlling the vor- 
tex structured At low temperatures, where the measure- 
ments are carried out, however, the distinct length scales 
likely survive, although they arc modified by the strength 
of the interband coupling, see below. Possible anisotropy 
of the gap on one or more Fermi surface sheets compli- 
cates the picture even further. We show here that judi- 
cious use of the quasiclassical approximation even with 
simplifying assumptions about the vortex structure can 
provide a general framework for the description of this 
problem, and a semiquantitative understanding of the 
new data on the BaFe2(Asi_ x P x )2 system. We compare 
our results with those obtained by Doppler-shift meth- 
ods, and show that if properly implemented this method 
also gives reasonable qualitative results in the low field 
range. 

This joint theory-experiment paper is organized as fol- 
lows. We first present our experimental results on the 
BaFe2(Asi_ :c P a: )2 system in Section HTl and compare to 
our previous results, as well as to data by other groups 
on the related heterovalcntly doped Ba-122 materials. In 
Section Unl we discuss the two-band quasiclassical model 
we use to study the system, and in Section llVl we give our 
results. Finally in Section [V] we present our conclusions. 



II. EXPERIMENTAL 

Tiny platelet crystals of BaFe2(Aso.7Prj.3)2 were pre- 
pared as described in Ref. 0- Subsequent measurements 
on crystals extracted from various positions in the cru- 
cible using x-ray diffraction and energy dispersion (EDX) 
analysis give a phosphorous concentration of 32.9±0.4%. 
A further test of the homogeneity of the crystals from a 
given growth batch is the measurement of the susceptibil- 
ity at the superconducting transition as shown in Fig. I 
of Ref. [HI. Here, for a collage of ~ 150 mg of these crys- 
tals a rather narrow transition was observed. A collage of 
18 mg of these microcrystals was then mounted on a sap- 
phire disk using GE7031 varnish. Approximately 75% of 
the crystals had the magnetic field perpendicular to the 
a-b plane (the plane of the crystals), whereas the remain- 
ing crystals were randomly oriented. The sapphire disk 
was mounted in our time constant method calorimeter jii 
and the specific heat from 0.4 to 7 K in fields from to 
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FIG. 1. (Color online) The original specific heat data*^ on 
BaFe2(Aso.7Po.3)2 as a function of field up to 15 T (solid 
symbols) with data from present work between 15 T and 
35 T (open symbols). Note the agreement between the linear, 
C/T oc H, extrapolation of the 15 T (colored finest) and 
35 T (black lines, present work) results. We extract 7 from 
the data using two (equivalent) methods: (a) by making an 
extrapolation C/T = 7 + /3T 2 + ST 4 from 2 K and above or 
(b) by taking the value of C/T at 1.5 and 2 K. The tempera- 
ture restriction eliminates both the influence of the anomaly 
and the field- induced nuclear contribution (see text), negligi- 
ble for H < 4 T above 1 K. The absolute accuracy of these 
data is ±5% while the precision of the data is approximately 
±2%. In addition, additional data with finer gradations in the 
measured fields up to 4 T were taken to explore the low field 
non-linear behavior. These data are shown on an expanded 
scale in Fig. [2] 



35 T was measured. Additionally, the specific heat of 
a standard (high purity Au) was measured in fields up 
to 14 T. Results on the standard (not shown) indicate 
agreement with published values to within ±3% in all 
fields. 



A. Results and Discussion 

The specific heat coefficient 7 = C/T of 
BaFe 2 (As . 7 Po.3)2 for < H < 35 T is shown by 
the open triangles in Fig. [1] There is a small low tem- 
perature anomaly in the specific heat data below about 
1.4 K (discussed in detail in[l3|). Such anomalies have 
been observed in other FcPn samples^ and in some 
cases, e.g. in BaFe2- a: C0; C As2, they show a rather strong 
magnetic field dependence.— However, as discussed in 
our previous report^ of the data up to 15T, the anomaly 
in BaFe2(Aso.7Po.3)2 is approximately field independent. 
Note that the small anomaly in the specific heat appears 
to vanish above 1.4 K, i.e. does not affect the estimate 
for 7 shown in Figs. [T] and [2] using data from 1.5 K and 
above. 

In order to have a closer look at the low field depen- 
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FIG. 2. (Color online) Low field 7 data up to 4 T from Fig.Q] 
on an expanded scale for T—2 K (blue), 1.75 K (red) and 1.5 K 
(black symbols). Green symbols are asymptotic limT->o C/T 
determined over the range 1.5 K< T <5 K. The fitting func- 
tions of the data are labeled beside the curves. Best power 
law fits to field dependence are shown in each case. 



dence of the specific heat, these data are shown on an 
expanded scale in Fig. [2] In our analysis below, we focus 
on the asymptotic T — s- behavior since it is directly 
related to the density of states at the Fermi level, which 
is easy to calculate reliably, and since it gives essentially 
the same field dependence as the nonzero T data. 



III. MODEL 



A. Quasiclassical approximation 



The quasiclassical (Eilcnberger) approximatioi 
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a powerful tool to describe the electronic properties of 
the superconducting state on the scales large compared 
to the lattice spacing, provided the condition kp£, 3> 1 is 
satisfied. Here kp is the Fermi momentum and £ the co- 
herence length. Since in this limit we can think of quasi- 
particles as propagating coherently along a well-defined 
trajectory in real space, this method is particularly well 
suited to address the inhomogeneous situations, such as 
the vortex state of type-II superconductors (SCs). An 
alternative and frequently used approach to the vortex 
state is to take into account the (classical) shift of the 
quasiparticle energy due to the local supercurrent flow. 
Such an approximation, often referred to as the Dopplcr- 
shift approach, is valid for nodal SCs with considerable 
weight of extended quasiparticle excitations ouside the 
vortex cores. Using this method, Volovik showed that for 
superconductors with line nodes these extended quasi- 
particle excitations lead to a non-linear magnetic field 
dependence of the spatially averaged residual density of 
states N(uj = Q,H) tx N oy / H/ H c2 , the result known as 
the Volovik effect^ This behavior was first confirmed by 



measurements of the specific hea t 23 i 24 and by subsequent 
calculations within the quasiclassical approximation for 
both a single vortex in a d-wave S C 25 ' 26 and for a vortex 
latticej 27 i 28 Both quasiclassical and Doppler-shift meth- 
ods fail at the lowest temperatures due to quantum ef- 
fects^, but in known systems with T c <C Ep these effects 
are negligible in practice. Both methods have successfully 
explained at a semiquantitative level the magnetic field 
dependence of the specific heat and thermal conductivity 
in a wide variety of unconventional superconductors^ It 
was also shown that the accurately calculated quasipar- 
ticle excitation spectrum is consistent with STM studies 
of the electronic structure around a vortex core. 2 ' 

Many experimental techniques which are sensitive to 
the low-energy density of states, such as thermal con- 
ductivity, specific heat, and NMR relaxation rate, can be 
used to draw conclusions about the possible existence and 
the momentum dependence of quasiparticle excitation in 
the bulk of iron-based superconductors (FeSCs) and thus 
about the structure of the superconducting gap and the 
distribution of gap nodes. The low T limit of the Som- 
merfeld coefficient in an applied magnetic field, 7 (if), is 
directly proportional to the spatially averaged local den- 
sity of states (LDOS) at the Fermi level. The Doppler- 
shift method has been used to calculate the LDOS for 
a two-band SC with two isotropic gaps of unequal size 
As =/= and to give an interpretation of the experimen- 
tal data available at that timeJ^ However, the Doppler- 
shift approach cannot account properly for the contribu- 
tions from the states in the vortex core that have a very 
large weight in the net DOS and hence gives a quantita- 
tively and sometimes qualitatively inaccurate description 
of the electronic structure of the vortex. For example, in 
a simple d-wave superconductor the spatial tails of the 
low-energy density of states around the vortex are aligned 
in the wrong directions.— To obtain a quantitative fit to 
the specific heat data presented in the previous section 
and to allow for a more decisive conclusion about the gap 
structure of BaFe2(Asrj.7Po.3)2, we will therefore use the 
quasiclassical approximation, which wc will briefly review 
in the following paragraphs. 

In the quasiclassical method, the Gorkov Green's func- 
tions are integrated with respect to the quasiparticle en- 
ergy measured from the Fermi level. The normal and 
anomalous components g(r, iw n ) and f(r,0,iu) n ) of 
the resulting propagator g obey the coupled Eilcnberger 
equations 



(la) 



2 [iu n + -v F ■ A(r)J + ihvp ■ V /(r, 9, iuj n ) 
= 2ig(r,6,iu n )A(r,6), 



2 \iuj n + -v F ■ A(r)J - ihvp ■ V f(r,6,ioj n 
= 2ig(r,6,iu n )A*(r,6), 



(lb) 



that have to be complemented by the normalization con- 
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1. 



(2) 



Here A(r, 0) is the order parameter, A(r) the vector po- 
tential, Vf is the Fermi velocity at the location at the 
Fermi surface labeled by 9, and w„ = (2n + l^ksT 
are the fermionic Matsubara frequencies. For two- 
dimensional cylindrical Fermi surfaces such as considered 
below, Vf = vpii where k = (cos 9, sin 9) and 9 is the an- 
gle measured from the [100] direction. In that case it is 
natural to write the position vector in cylindrical coordi- 
nates, r = (p, <fi, z), where <j> is the winding angle around 
the vortex in real space. 

Making use of the symmetries^ of the quasiclassical 
propagator 



:i:s 



/(r,k F ,iu; n ) = f*(r,\t. F ,-iu n ), (3a) 
/(r, -k_p, -iui n ) = /(r, k F , ioj n ), (3b) 
g(r,k F ,iw n ) = g*(r, k F ,-ioj n ), (3c) 

the diagonal part of the normalization condition ([2]) can 
be written in a more explicit form as 

[g(r, 9, ico n )] 2 + /(r, 9, iw n )/* (r, 9 + tt, iw n ) = 1 . (4) 

Instead of solving the complicated coupled Eilenberger 
equations everywhere in space, we follow Refs. [26| and 
|32| and parameterize the quasiclassical propagator along 
real space trajectories r(x) = r + xv F by a set of scalar 
amplitudes o(x) and b(x), 



ff(r(a:)) 



1 



1 + a(x)b(x) 



1 



- a(x)b(x) 
2b(x) 



2a(x) 
-1 + a(x)b(x) 



(5) 



These amplitudes obey numerically stable Riccati equa- 
tions, 

v F d x a{x) + [2Cj n + A* (x)a(x)}a(x) - A{x) = 0, (6a) 
v F d x b{x) - [2G) n + A(x)b(x)]b(x) +A*(x)=0. (6b) 

For the single vortex problem the spatial dependence van- 
ishes far away from the vortex core, and hence we have 
the initial conditions 



a(— 00) = 
b(+oo) = 



A(-oo) 



u n + v/w2 + |A(-oo)| 2 ' 

A* (+00) 
w„ + s/ujI + |A(+oo)| 2 ' 



(7a) 
(7b) 



Here we have set h = 1 and we have introduced the modi- 
fied Matsubara frequencies iui n (x) = iuj n + (e/c)v F -A(x). 
Since the modification of the Matsubara frequencies due 
to the external field is of the order of 1/k 2 where k = 
Xl/S, is the ratio of the London penetration depth and 
the coherence length the term proportional to A(x) in 



Eq. ([6]) can be neglected for strong typc-II superconduc- 
tors. 

After an analytic continuation of the Matsubara fre- 
quencies to the real axis, iuj n — > uj + iS, the local density 
of states can be calculated as the Fermi surface average 
of the quasiclassical propagator 
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ab 



(8) 



where No is the normal density of states at the Fermi en- 
ergy. To obtain stable numerical solutions we use a small 
imaginary part S = 0.02T C in the analytical continuation, 
where T c is the critical temperature of the superconduc- 
tor. 



B. Two-band model 

The Fermi surface of the optimally doped 
BaFe2(Aso.7Po.3)2 consists of multiple Fermi sur- 
face sheets. DFT calculations showed that there are 
three concentric hole cylinders in the center of the 
Brillouin zone (r point) and two electron pockets at the 
zone corner (X point) M- Laser ARPES measurements^ 
found a superconducting order parameter that is fully 
gapped with comparably sized gaps on each hole pocket 
of the order of A^/fegTc ~ 1.7. Taking into account 
the results from thermal conductivit y 10 ' 12 and NMR 
measurementsii as well as the measurements of the 
specific heat coefficient in low fields presented above, 
that all consistently report evidence for low-energy 
quasiparticles, this ARPES result implies a nodal gap 
on the electron pockets. 

For numerical convenience we adopt below a two-band 
model, distinguishing only between electron and hole 
pockets. Inclusion of all Fermi surface sheets then only 
enters as a weighting factor for the electron and hole 
pocket contributions as we discuss in the following sec- 
tion. We take the gaps on the electron and hole pockets 
in the form Ai,2(0) = Ag' ?l <I>i i 2(6'), where the angle 9 pa- 
rameterizes the appropriate Fermi surface, assumed to be 
cylindrical. We assume an anisotropic gap on the elec- 
tron pocket $i(0) = (1 + rcos26»)/ N /l + r 2 /2, and an 
isotropic gap around the hole Fermi surface, $2(0) = 1- 
If the anisotropy factor r > 1, the superconducting gap in 
the electron band, Ai(0), has accidental nodes; if r = 0, 
Ai(9) is isotropic like Az^)- 

First we assume Aq — Aq, as is often found by ARPES 
(at this writing there are no ARPES results on this mate- 
rial which resolve the gap on the electron pocket). Since 
we consider well separated electron and hole bands, we 
can solve the Riccati equations, Eqs. ([6]), for the two 
propagators separately, and the only coupling of the 
pockets is via the self-consistency equations on the or- 
der parameter, see below. With this in mind we nor- 
malize the energy and length for the electron and hole 
bands by the gap amplitudes Aq and Aq, and the co- 
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FIG. 3. (Color online) The spatially averaged ZDOS, normal- 
ized to the normal state value N(u) — 0)/No for a nodeless 
(orange) and a nodal (blue) single-band superconductor. The 
dashed lines show the idealized linear H and y/~H behavior 
for a clean s-wave and d-wave SC, respectively. The symbols 
are numerical results for a single band SC with an isotropic 
s-wave gap (circles) and a strongly anisotropic nodal gap (tri- 
angles). Additionally we compare results with (solid symbols) 
and without (open symbols) taking into account the vortex 
core reduction due to the Kramer-Pesch effect. Here we have 
ignored the field dependence of the superconducting gap, i.e. 
A(H) = A . 

herence lengths £q = v f /Aq and £q = v f /Aq respec- 
tively. Fermi velocities therefore appear as an input. 
DFT calculations for a comparable Ba-122 system 3 -^ give 
v h F = 1.979 x 10 5 m/s and vp = 3.023 x 10 5 m/s, i.e. 
vp/vp = £q/£o = 0.65. In our analysis we keep this 
ratio but reduce the value of both Fermi velocities by a 
factor of 5 to approximately account for the mass renor- 
malization of this system near optimal doping!*^ This 
reduction also gives a roughly correct value of the c-axis 
upper critical field H C 2 ~ 50 T. In the limit of negligible 
coupling between the bands, the upper critical field is 
determined by the overlap of the vortices with smallest 
core size, 

* = JL = A /?£ (9) 

Below we solve the Eilenberger equations and determine 
the density of states for an isolated vortex and for each 
band separately In a two-band system the spatial pro- 
file of the quasiparticle states on the electron and hole 
bands is controlled by the respective coherence lengths, 
and therefore spatial averaging weighs the contributions 
of the bands differently compared to the DOS of a sys- 
tem with a single or two equal coherence lengths. This is 
the most significant difference compared to a single-band 
model. 

The superconducting order parameters in the two 
bands are related by the interband component of the 



pairing interaction. We consider a general coupling ma- 
trix in the factorized form, \ vll {9,9') = \ VI1 $ V (#)<I> M (#'), 
where p, v = 1,2 and = V^^N^. Here Vn = V e and 
V22 = Vh are the intraband pairing interactions in the 
electron and the hole band, respectively, while V12 = V e % 
is the interband interaction. is the normal density of 
states at the Fermi level. Then the gap equation for an 
inhomogeneous superconductor is 

A u (r) = 2irT £ £ (M0)/ M (r, 6, iu n ))g . (10) 

/i=l,2 lo„>0 

Here A„(r) is the momentum independent part of the 
gap function; Ai )2 = A e ' h at T = and H = 0. 

In the vortex state the self consistent determination 
of the spatially dependent order parameter is a complex 
task. Since we are interested in relatively low fields, when 
the vortices are well separated, we solve the Eilenberger 
equations for the order parameter that is assumed to have 
a single vortex form, 

A-(f H;0) = A^mtanh ( 1 + rcos20 (lla) 

A h (p,H) = A 2 (i/)tanh (^^j . (lib) 

Here p = (p, 4>) is the two-dimensional projection of the 
radius vector in cylindrical coordinates, and factor of 0.1 
is introduced to approximate the shrinking of the core 
size in the self-consistent treatment at low temperatures 
(Kramer-Pesch effec t 39 ' 40 ) . This single vortex ansatz 
provides a qualitatively correct description of the low- 
field state, close to what is found by full numerical solu- 
tion^ To account for the suppression of the bulk order 
parameter by the magnetic field, we determine the coef- 
ficients Ai^{H) from the Pesch approximation^ where 
in the presence of an Abrikosov lattice the diagonal com- 
ponents of the Green's function by its value averaged 
over a unit cell of the vortex lattice. This approximation 
proven to give reliable results over a considerable range 
of magnetic fields and is incorporated into our approach. 

Note that our ansatz for the order parameter becomes 
quantitatively inaccurate for strong interband coupling in 
the regime of applicability of the Ginzburg-Landau the- 
ory since the core sizes of the two bands approach each 
other42 We verified in a fully self-consistent calculation 
that in the parameter range that we use the correspond- 
ing effect on the specific heat is of order 1% or less and 
hence can be neglected. We therefore use Eq. (fl~Tj) here- 
after. 

To proceed we substitute Eq. (fTTj) into Eq. ©, solve 
for a(x) and b(x), and use Eq. © to find the local den- 
sity of states N(p,H). To approximate the specific heat 
coefficient, we evaluate the spatial average of the zero 
energy local density of states 

N(H) = d<t> [ R dppWp-, (12) 
Jo Jo ttR 2 N 



6 



TABLE I. The different models for the coupling matrix and 
the gap anisotropy on the electron pockets considered in this 
work. 



An 



Ai 



A21 



A 2 



T c /K H c2 /T 



case (1) 0.51 0.51 0.33 0.65 0.9 31 54 

case (2) 1.00 0.02 0.013 0.81 0.9 31 47 

case (3) 0.51 0.51 0.34 0.64 1.3 31 54 

case (4) 1.00 0.023 0.015 0.77 1.3 31 42 



where the intervortex distance R depends on H as de- 
scribed by Eq. ©. The total density of states is then 



given as 



N{H) t 



w e N e (H) +w h N h (H) 



Wh 



(13) 



where w e /wh = 2N§/Nq = 2£ if we consider, for exam- 
ple, two electron Fermi surface sheets in the folded Bril- 
louin zone and denote ( = Nq/Nq = v F /vp = 0.65 = 
A21/A12. 



IV. RESULTS 

To illustrate that the salient features of the vortex state 
DOS are captured in our approach in Fig. [3] we show the 
field dependence of the spatially averaged zero energy 
local density of states (ZDOS) for a one-band SC with 
either an isotropic s-wave gap or a strongly anisotropic 
nodal gap (r = 1.3). Note that, while the field depen- 
dences in both the nodal and fully gapped cases clearly 
fit the anticipated power laws at low fields, \Jli and if, 
respectively, there is a significant influence on the magni- 
tude of the DOS caused by the size of the core, with the 
smaller core size yielding smaller ZDOS. In particular, 
in the absence of the Kramer-Pesch effect, for the nodal 
case the ZDOS would exceed the normal state value at 
fields far below H C 2, which is unphysical. 

Below we consider r = 0.9 and r = 1.3 to mimic a gap 
with deep minima and accidental nodes, respectively. To 
show different types of behavior allowed within our mi- 
croscopic model we chose four sets of coupling constants, 
two for each value of r, as shown in Table HI In cases (1) 
and (3), the interband pairing A12 is strong and close to 
the intraband parameter An, while in case (2) and (4) 
A12 <C An, A22- 

In Fig. 2ta) we show the self-consistently determined 
magnitudes of the bulk gaps in the vortex state Ai^-ff) 
as defined in Eq. (TO]) and (JTTJl- H c2 ~ 40 - 50 T. 
In the cases with only weak interband pairing (2) and 
(4), the gap on the electron Fermi surface deviates 
considerably from the phenomcnological form A (if) = 
Aoa/1 — H/H C 2- Figs. SJb) and (c) show the spatially 
averaged ZDOS corresponding to each band. For N e {H) 
and for r = 1.3 the \/H behavior of the Volovik effect 



is clearly visible at lower fields. Comparing Fig. @|b) 
to Fig. [3] we find that within the two-band model the 
density of states on the electron band N e (H) reaches a 
quasi-linear behavior already at smaller fields than the 
corresponding density of states for the one-band case. In 
Fig.[3]a linear behavior is never observed, and might only 
be fit over some intermediate field range for H/H C 2 > 0.2, 
while in the multiband case N e (H) displays a clear linear 
behavior already for H/H C 2 > 0.1. 

It is tempting to interpret the low-field crossover to a 
quasilinear field variation as evidence for a small energy 
scale A sm = Ag(l — r)/yl + r 2 /2 on the electron band; 
this, however, seems unlikely. Provided A sm <C Ag, the 
gap still increases linearly along the Fermi surface away 
from the nodal points above this energy scale, simply 
with a different slope. Then within the usual Volovik 
argumentation the contributions from extended states at 
these intermediate energies give rise to a s/H contribu- 
tion even if A sm < Eh <C A maa; , where Eh oc y/~H is the 
average Doppler shift and A max = Ag(l + r)/y/l + r 2 /2 
is the maximum gap. There is therefore no true linear- H 
behavior arising from the electron band with gap nodes. 
Consequently, we interpret this crossover as the conse- 
quence of the two-band behavior coupled with a gradu- 
ally increasing contribution of core states which is nearly 
linear in field. Fig. |4jc) clearly shows that the density of 
states on the hole band N h (H), assumed here to be fully 
gapped, is always linear as a function of field and the 
results for the two different coupling matrices considered 
here are very similar. However, as mentioned before, the 
slope is smaller than the one predicted for an idealized 
s-wave SC. 

Using Eq. (fT3"|) . the spatially averaged ZDOS on the 
electron and the hole band are added with different 
weights. Using the results presented in Figs. @|b) and 
(c) as case (4), we investigate several cases. Since there 
are two electron pockets, and assuming that only one hole 
pocket contributes significantly to the low energy density 
of states (or that a naive average over the hole pockets is 
sufficient), the net DOS and the field dependence of the 
Sommcrfeld coefficient are only functions of the ratio of 
the densities of states on the electron and hole sheets. In 
the following we will study three cases that we will ab- 
breviate with "Q" indicating the use of the quasiclassical, 
or Eilcnbcrgcr, approach: 

• (Qa): we assume that only one hole pocket con- 
tributes considerably to the low energy DOS, and 
use the weights w e /wh = 2Nq/Nq taken from the 
DFT calculation, N§/Nft = 0.65, see Ref. US 

• (Qb): We once again fix Nq/Nq = 0.65, but adopt 
a model for which the normal DOS for all three hole 
pockets of Ba2Fc2(Aso.7Po.3)2 are the same and for 
which all three pockets contribute equally to the 
low energy DOS, hence w e /w h = 2N^j/3N^; 

• (Qc): Wc do not hold the ratio N§/N^ fixed, 
but instead calculate the weights for the electron 
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FIG. 4. (Color online) Results of quasiclassical calculations for the parameters in Table I. (a) Magnetic field dependence of the 
gaps in the two-band model calculated within the Pesch approximatio n 41 ' 43 for case (1) - (4). We assume A e (H = 0) = A h (H = 
0) here. The four sets of coupling constants Xij are listed in TABLE U (b) Field depend ence of the space average ZDOS N e (H) 
on the electron pocket for the four cases with anisotropic gap with angular variation $ e (0) = (1 + r cos 20)/ y/ 1 + r 2 /2. (c) 
Field dependence of the space average ZDOS N h (H) for the four cases with isotropic gap along the hole pocket. 




FIG. 5. (Color online) Comparison of the experimentally mea- 
sured normalized specific heat coefficient (large pink dots) to 
different theoretical results for the spatially averaged ZDOS. 
The dotted violet and solid orange curves are the predic- 
tions for the spatially averaged ZDOS for a clean s-wave 
and d-wave SC. The blue squares (case (Qa)) and green 
diamonds (case (Qb)) are the differently weighted sums of 
N e (H) and N h {H) evaluated for case (4) of Figs, g^b) and 
(c). The black line (case (Qc)) is obtained using the formula 
7 tot = aiN e (H) + a 2 N h (H) where at = 3.2 mJ/(mole K 2 ), 
a? = 10.3 mJ/(mole K 2 ) are determined with the least square 
fit to experimental data below 30 T. Note "d-wave" and "s- 
wave" curves represent simple extrapolations of the low-field 
s/H and H terms up to H C 2- 



pockets a\ and for the hole pockets 02 by a least 
squares fit to the experimental data using the for- 
mula 7 tot = aiN e (H) + a 2 N h (H). If we normal- 
ize it to the presumed contribution of the super- 
conducting fraction, j n — 70 ~ 14 mJ/(mole K 2 ), 
where 70 is the extraneous term, see below, we find 



w e /(w e + w h ) = ai/(7n ~ 7o ) and w h /(w e + w h ) = 
02/(7™ - 7o) and a 1 /a 2 = w e /w h . 

In Fig. [S] wc compare the results for all three cases to the 
experimentally measured specific heat coefficient (pink 
dots). The experimental values are obtained by extrap- 
olating the measured specific heat coefficient 7 at var- 
ious temperatures to T = 0. The upper critical field 
H C 2 is taken to be 52 T, see Ref. [l2|. The normal state 
7,, = 16 mJ/(mol K 2 ) can be obtained by extrapolating 
7 to H C 2- A substantial residual^ 70 = 1.7 mJ/(mol K 2 ) 
in the superconducting state, presumed due to disorder, 
is subtracted in the plots of the field dependence from 
the experimental data (pink dots) to compare with our 
quasiclassical calculation in the clean limit (blue squares 
and green diamonds). Note that subtracting of the resid- 
ual C/T tends to enhance the scatter in the low-T data 
of Fig 2. 

From Fig. [SJ we see that the results derived for model 
(Qb) with three equal mass hole pockets and two equal 
mass electron pockets are in good agreement with the 
experimental data: both experiment and theory show a 
"Volovik effect" at the lowest fields and then a crossover 
to a linear H dependence at intermediate fields. While 
model (Qa) has the same qualitative behavior, the rel- 
ative weights of hole and electron bands are apparently 
not consistent with the normalized experimental data, 
and the fit is much poorer. Compared to (Qb) the least 
squares fit (Qc) to the experimental data (black line) is 
only marginally improved, and gives Nq/Nq = 0.47 with 
two electron pockets/three hole pockets or 0.16 with two 
electron pockets/one hole pocket, same order as obtained 
from DFT. 

For completeness it is important to determine whether 
the experimental data can be appropriately fit within the 
confines of a simple two-band Doppler shift approach. 
We detail this method in the Appendix, where we show 
that models (a) and (b) do not give a satisfactory fit to 
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FIG. 6. (Color online) 7 tot plotted with experimental data in 
absolute units. Case (Qc) is the similar black line in Fig. [S] 
(Dc) is obtained using formula j tot = aiN e (H) + a2N h (H), 
where N e (H) and Nh{H) are represented by the open squares 
and circles in Fig. [8] and ai = 1.50 mJ/(mole K 2 ), ai = 
65.6 mJ/(mole K 2 ) for (Dc) are determined with the least 
square fit to the experimental data below 15 T. 
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FIG. 7. (Color online) Experimentally measured specific heat 
coefficient (large pink dots) compared to calculations with 
deep gap minima (case (1), r = 0.9, blue triangles) and acci- 
dental nodes (case (4), r = 1.3, green squares). In both cases 
the weight of electron and hole pocket contributions has been 
chosen in agreement with case (Qb). 



the experiment. In contrast, model (c) yields a rather 
similar field dependence of the field-induced enhance- 
ment of the Sommerfeld coefficient for the quasiclassi- 
cal and Doppler (Dc) methods, as shown in Fig. [5] At 
the same time the best fit linear coefficients for (Dc), 
ai = 1.50 mJ/(mole K 2 ), a 2 = 65.6 mJ/(mole K 2 ), give 
the ratio of the normal state DOS for two electron/three 
hole Fermi sheets of N^/Nf} = Sa 1 /2a 2 « 0.03, very dif- 
ferent from the value of 0.65 obtained within the band 
structure calculations. Consequently, the quasiclassical 
methods provides a far more satisfying fit to the data. 

As is usually the case with the measurements that 
probe the amplitude rather than the phase of the gap, it 
is difficult to distinguish the deep minima from the nodes. 
In this case we find that with our current uncertainty in 
the band parameters, and the scatter in the data, it is 
impossible to assert the nodal behavior purely from the 
current data. Fig. [7] shows the comparison of cases (1) 
and case (4) of Table U corresponding to r = 1.3 and 
0.9, i.e. with and without true nodes, with the weights 
of case (Qb). Even though the nodal fit appears bet- 
ter at the lowest fields, higher H data are in between 
the two cases. Therefore the conclusion about the true 
node comes from the data on other experiments, such as 
penetration depth. 



V. CONCLUSIONS 

Among the various families of Fe-based superconduc- 
tors, BaFc2(Asi_ x P 2; )2 may be a key system for under- 
standing the origins of superconductivity. In part this is 
because, alone among the materials thought to display 



nodes in the superconducting gap, it possesses a rather 
high T c of 31 K, and hence the interplay of the pairing 
mechanism and Fermi surface shape and parameters in 
determining the gap anisotropy is under special scrutiny. 
The lack of an observable Volovik effect in earlier spe- 
cific heat measurements was a cautionary note in an oth- 
erwise consistent array of measurements in support of 
gap nodes. In this paper, we have presented new exper- 
imental data at both lower and higher fields than previ- 
ous measurements, and found that the initially reported 
linear-lf behavior extends up to 35 T, but that at low 
fields (H < 4 T) more precise measurements with smaller 
gradations in the change of field between data points are 
now clearly consistent with a Volovik-typc effect. The 
residual T — > Sommerfeld coefficient ^(T — > 0) is 
about 1.7 mJ/mol-K 2 , consistent with possible nanoscale 
disorder in the sample. The low-field sublinear depen- 
dence of the Sommerfeld coefficient is a strong indication 
that nodes (or deep minima) are present, and provides 
the sought-after consistency with other probes without 
having to make extreme assumptions about the ratio of 
masses on electron pockets to those on hole pockets, as 
was proposed in Ref. [H|. 

It is nevertheless striking that indications of nodal be- 
havior on the same samples is so much weaker in the 
specific heat measurements as compared to thermal con- 
ductivity and penetration depth. This is clearly indi- 
cating that the nodes are located on the pockets with 
smaller masses and/or longer lifetimes, as was pointed 
out in Ref. [l3|. We have attempted to put this statement 
on a semiquantitative basis by presenting a quasiclassi- 
cal (Eilcnbcrger) calculation of the density of states and 
specific heat of a two-band anisotropic s± superconduc- 
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tor. Comparison with the Dopplcr shift method allowed 
us to argue that the quasiclassical calculation is superior 
for semiquantitative purposes. We find that the unusu- 
ally small range of Volovik-type behavior, followed by a 
large range of linear- H behavior, is due to the small gap 
and weak nodes on the small mass (presumably electron) 
sheeti 12 ' 13 Good fits to the data are obtained for aver- 
age hole and electron maximum gaps of approximately 
equal magnitude, in the weak intcrband coupling limit. 
The success of this fit should not, however, tempt one 
to draw definitive conclusions about the relative magni- 
tudes of the pairing interactions. The proliferation of 
parameters in the theory make it difficult to determine 
gap magnitudes, density of states ratios, and nodal prop- 
erties with any quantitative certainty. Equally good fits 
can be obtained, for example, with substantially smaller 
full gaps than anisotropic gaps; the nodes control the low- 
field behavior, and the small full gap gives rise to a large 
linear term. What is important is that we have shown 
that a fit can be obtained, with reasonable values of the 
parameters, that it can only be obtained if nodes exist 
on one of the Fermi sheets, and that it requires going be- 
yond the simple Doppler shift picture. It is our hope that 
the results of this calculation and fit will eventually lead 
to a more quantitative first principles based calculation. 
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Appendix A: Comparison with the Doppler-shift 
method 



In the following we briefly discuss the basic concepts 
of the Doppler-shift method and compare it to the quasi- 
classical approximation as manifested in the formulation 
of the Eilenbcrgcr equations introduced in the main text. 
The Doppler shifted energy due to the local supercurrent 
flow is lj — mvp ■ v s (r) where 



e,h 



v.(r) = (0) 
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2|r| 



■ sin 
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2p 



ne,h 
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Here, we assume A e (T = 0,H = 0) = A h (T = 0,H = 
0) = Aq h and use p e ' h = |r|/£o' • Therefore the normal- 



ized local DOS is 
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and thus the normalized spatially averaged DOS is 
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(A3) 



Here we have introduced the normalized vortex cell ra- 
dius R = R/£q for the electron bands or i?/£o for the 
hole bands, respectively. Since the Doppler-shift method 
does not capture core state contributions it underesti- 
mates the slope of the magnetic field dependence of the 
zero energy DOS of an s-wave SC. Since the core region 
only gives negligible contributions to the total DOS one 
can in principle avoid the divergence of the Doppler-shift 
energy as p — > by cutting out the complete core region 
with a lower limit £o for the radial integration. Here we 
have included the core region when integrating over the 
vortex unit cell. To model A(9) we use a similar function 
as given by Eq. but without explicitcly modeling the 
core structure 



A e (ft6)=A 1 (H = 0) 
A h (p)=A 2 (H = Q), 



1 + r cos 26 
\/l + r 2 /2' 



(A4a) 
(A4b) 



and we use the self-consistently calculated Ai^-ff = 0) 
in case (4) of the quasiclassical calculation in which the 
anisotropy factor r = 1.3 for the gap along the electron 
Fermi surface sheet and the ratio of the normal DOS at 
the Fermi energy is taken as N§/Nq — 0.65. In Fig. [5] 
we show the results obtained within the Doppler-shift 
approach and compare them to case (4) of the quasi- 
classical method. Again we also show predictions for an 
idealized clean s- and d-wave SC. We conclude that the 
Doppler-shift method and the quasiclassical method give 
comparable results at lowest fields but start to deviate 
as soon as the field increases. One reason might be that 
with increasing field and decreasing inter- vortex distance 
the core states that are not correctly accounted for within 
the Doppler-shift method but captured within the quasi- 
classical approach become increasingly more important. 
However, due to the limitations of the single vortex ap- 
proximation the overlapping of vortices is not correctly 
reproduced and the DOS is overestimated (In Fig. [5] the 
blue triangles rise too fast). 

In Fig. [5] we compare the least squares fit by Doppler- 
shift method (blue dashed curve) together with least 
square fit by quasiclassical method (the similar black 
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FIG. 8. (Color online) Normalized specific heat coefficient 
determined within the quasiclassical approach (blue triangles 
and blue dashed curve) and within the Doppler-shift method 
(open squares and circles). Total normalized specific heat 
coefficients obtained within Doppler-shift method are shown 
by magenta dots, olive dash-dots and black solid line for cases 
(Da, Db, Dc), corresponding to (Qa, Qb, Qc), respectively. 
The dotted violet and solid orange curves are the power laws 
predicted for the spatially averaged ZDOS of an idealized s- 
wave and d-wave SC. Large pink dots show the experimentally 
determined normalized specific heat coefficient. 



line from Fig. [5]) and experimental data (large pink 
dots). The linear coefficients for (Dc) are a\ = 
1.50 mJ/(mole K 2 ), 0,2 = 65.6 mJ/(molc K 2 ). Compared 
to the linear coefficients for (Qc) (ai = 3.2 mJ/(mole K 2 ), 
a-2 = 10.3 mJ/(mole K 2 )), they give a nonphysical ratio 
of the normal DOS at Fermi energy if we consider two 
electron Fermi sheets and three hole Fermi sheets. To see 
this point, let's consider equation 



C e , h (T,H) 



+ ce- 



ded 



LU 2 N e , h (H,L0) 



T 2 cosh(w/2T) 
AN e<h (H,u)T (asT^O) 



(A5) 



7e ,h « AN e , h (H, 0) = AN°' h N e , h (H, 0). 

Here A is a numeric constant and we write N et h(H, 0) = 
No' h N e!h (H,0) where N e , h (H,0) is the ZDOS calculated 



by the Green's function method as defined in Eq. (|12[) . 
Denote the number of Fermi sheets included in summa- 
tion as rips and define npgANfj = ai and ii fs ANq 



1 fs N o_ 



and Wh 



FS 



N$). 



OL2 (equivalent to w e = n'j 

Therefore 7 tot = a x N e {H) + a 2 N h (H). Note N e . h {H) 
are dimcnsionless and are in unit of mJ/(molcK 2 ). 
Optimized parameters ai t 2 for least square fit (Dc) to 
the experimental data below 15 T are a\ = 1.50 and 
«2 = 65.6. This leads to our estimate in the main text 
of oi/aa = (n FS ANS)/(n FS ANft) = (2NS)/(3N>}), and 
Nq/Nq = 0.03, and to our conclusion that Doppler-shift 
method does not provide a satisfying physical explana- 
tion to our specific heat experiment. 
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